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ABSTRACT 

Adaptive spatial domains are currently used in Smooth Particle Hydrodynamics (SPH) 
with the aim of performing better spatial interpolations, mainly for expanding or shock 
gas dynamics. In this work, we propose a SPH interpolating Kernel reformulation suit- 
able also to treat free edge boundaries in the computational domain. Application to 
both inviscid and viscous stationary low compressibility accretion disc models in Close 
Binaries (CB) are shown. The investigation carried out in this paper is a consequence 
of the fact that a low compressibility modelling is crucial to check numerical reliability. 
Results show that physical viscosity supports a well-bound accretion disc formation, 
despite the low gas compressibility, when a Gaussian-derived Kernel (from the Error 
Function) is assumed, in extended particle range - whose Half Width at Half Max- 
imum (HWHM) is fixed to a constant h value - without any spatial restrictions on 
its radial interaction (hereinafter GASPHER) . At the same time, G ASPHER ensures 
adequate particle interpolations at the boundary free edges. Both SPH and adaptive 
SPH (hereinafter ASPH) methods lack accuracy if there are not constraints on the 
boundary conditions, in particular at the edge of the particle envelope: Free Edge 
(FE) conditions. In SPH, an inefficient particle interpolation involves a few neighbour 
particles; instead, in the second case, non-physical effects involve both the boundary 
layer particles themselves and the radial transport. 

Either in a regime where FE conditions involve the computational domain, or in a 
viscous fluid dynamics, or both, a GASPHER scheme can be rightly adopted in such 
troublesome physical regimes. 

Despite the applied low compressibility condition, viscous GASPHER model shows 
clear spiral pattern profiles demonstrating the better quality of results compared to 
SPH viscous ones. Moreover a successful comparison of results concerning GASPHER 
ID inviscid shock tube with analytical solution is also reported. 

Key words: accretion, accretion discs ~ hydrodynamics - methods: numerical, N- 
body simulations - stars: binaries: close, dwarf novae, cataclysmic variables 



1 INTRODUCTION 



In its original version (|Monaghanl Il985l , 1 19921 : 
iMonaghan fc Lattanzid [l985l ') SPH was formulated adopt- 
ing a constant particle smoothing length (spatial smoothing 
resolution length or resolving power) h, where the adopted 
interpolation Kernel works, to perform free Lagrangian 
gas dynamics. ASPH methods are currently adopted with 
the aim of performing better spatial interpo lations mainly 
in expanding or in shock gas dynamics (jEvrardl 1 19881 : 
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In order to build up a well-bound accretion disc in invis- 
cid conditions, the ejection rate at the disc's outer edge must 
be at least two or three times smaller than the accretion rate 
at the disc's inner edge. Whenever this condition is fulfilled. 
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the disc's outer edge, as well as the whole disc, does not dis- 
perse in spite of high pressure forces which are also depen- 
dent on the gas compressibility: — Vp/p = —(7— l)V(pe)/p. 
Therefore, low compressibility gases are naturally more eas- 
ily sensitive to the loss of blobs of gas at the disc's outer edge 
itself, towards the empty external space, if the gravitational 
field is not able to keep disc gas in the gravitational po- 
tential well. Such effe cts are enhanced and stro ngly evident 
in inv iscid conditions l|Molteni et al.|[T99ll : iLanzafame et al., 
Il992l ) and the moderate contribution of artificial viscosity 
terms does not prevent such effects. Such a viscosity does 
not work like a true physical one since it operates only when 
different fiuid components approach each other, being zero 
during fluid particle repulsion. 

High compressibility gas dynamics does not allow us 
to distinguish the truth regarding whether a technique 
is able to perform a correct fiuid dynamics. In fact, in 
such a modelling, accretion discs woul d be formed anyway 
even in physically inv iscid conditions (|Molteni et al.|[l99ll : 
ILanzafame et aLlll992l ) and the role of Kernel choice and of 
its resolving power are hidden. To stress such an idea, in 
this work a low compressibility 7 = 5/3 polytropic index is 
adopted throughout, working with the same binary system 
parameters such as stellar masses and their separation and 
adopting the largest value {ass = l)as for the Shakura and 
Sunyaev viscosity prescription. 

In this paper, physically inviscid and viscid disc mod- 
els are shown, where a more suitable Gaussian-derived Ker- 
nel formulation, as far as both transport mechanisms and 
expanding or collapsing gas dynamics are concerned, is 
adopted. Throughout the accretion disc models, the same 
supersonic mass transfer condition at LI are adopted. 

The numerical scheme here adopted, as any other nu- 
merical method, is characterized by the assumed spatial 
smoothing resolution length h. The mass and angular mo- 
mentum radial transport is also affected by the SPH particle 
smoothing resolution length h. Too small h values prevent 
the radial transport, while large h values produce a too ef- 
fective radial transport of matter towards the centre of the 
gravitational potential well, as well as of angular momen- 
tum toward the disc's outer edge. A large h ensures a high 
particle overlapping (interpolation) but at the same time 
it produces a strong particle repulsion rate due to pressure 
forces especially in low compressibility regimes on the disc's 
outer FE. On the contrary, a too small smoothing reso- 
lution length h compromises any fluid dynamic behaviour 
and shock handling. The artiflcial viscosity term prevents 
spurious heating and handles shocks as a "shock capturing 
method" . The artiflcial viscosity is a function of the smooth- 
ing resolution length itself or of some kind of spatial length. 
A too small h value does not prevent particle interpenetra- 
tion, destroying any fluid behaviour, because of lack of artifl- 
cial viscosity,_Moltcru ct al. ( 1991); Lanzaf ame ct al. (1992); 
ILanzafame! (|2003l . 120091 ) and lLanzafame et al.l (|200^) discuss 
what we statistically define as a well-defined and bound ac- 
cretion disc. As far as the numerical resolution is concerned, 
a number of disordered neighbour particles of the order of 
10 (more or less) is considered, in principle, the minimum 
number of neighbours in order to achieve an adequate 3D 
numerical interpolation, although a number of neighbours 
larger than 30 is currently adopted to achieve a higher accu- 
racy. This is the criterion we adopted to define a well-bound 



accretion disc. Lesser neighbours for each particle are con- 
sidered an unsuitable number as far as both interpolation 
efficiency and disc binding into the primary's gravitational 
potential well are concerned. 

In the next sections, after discerning the artificial and 
the turbulent physical viscosities, we describe how ASPH 
techniques work and their limits when the viscous transport 
and/or FE conditions are involved, as well as why GAS- 
PHER could be a solution. In particular, in §2 we compare 
how artificial and turbulent physical viscosities differently 
work; in §3 we show how GASPHER works and why it 
does not suffer of some SPH and/or ASPH lack. At last, 
in §5 we report 3D accretion disc results showing some in- 
teresting features in our viscous simulations, whilst in §6 
we discuss on the accuracy of SPH-derived techniques. In 
the Appendix, after showing the mathematical background 
underlying SPH-derived schemes (for readers knowing how 
SPH and ASPH work, this mathematical section can be eas- 
ily skipped without any difficulty, being instead essential for 
others), we also compare results of GASPHER, SPH and 
ASPH non viscous ID and 2D selected tests. A comparison 
with analytical solutions is also given, whenever it is possi- 
ble. 



2 THE ARTIFICIAL AND THE TURBULENT 
PHYSICAL VISCOSITIES 

In our physically viscous disc modellin g, the Sh akura 
and Sunyaev prescrip tion (|Shakural Il972l . Il973l : 
IShakura fc SunvaevI Il973l ) is adopted with the largest 
ass = 1 value t o str e ss numerical reli ability of results 
(ILanzafame et al.l I2OO6I : iLanzafamd 120091 ). The SPH for- 
mulation of viscous contributions in the Navier-Stokes 
and en ergy equations has been developed by iFlebbe et al.l 
(|l994al lbl). These goals are not obtained by artificial viscos- 
ity which is, however, introduced in both models to resolve 
shocks numerically and to avoid spurious heating. Artificial 
viscosity vanishes when the limit value of the particl e 
interpola t ion d o main goe s to z ero. iMeglicki et al.l l ll993h : 
iDrimmell l| 19961 ): iMurravl (|l996l ) and lOkazaki et all (|200a ) 
demonstrated that the linear component of the artificial 
viscosity itself, in the continuum limit, yields a viscous shear 
force. In particular, the last two authors have explicitly 
formulated such an artificial viscosity contribution in the 
momentu m and ene rgy equations. Moreover. iMurravl (|l996l ) 
and Okaz aki et al.l (2002) found an analogy between the 
shear viscosity generated by the linear artificial viscosity 
term and the well-known Shakura and Sunyaev shear 
viscosity, in the continuum limit. SPH method, like other 
finite difference schemes, is far from the continuum limit; 
moreover we need the quadratic [Psph, Von Neumann- 
Richtmyer-like viscosity) artificial viscosity term to handle 
strong shocks. Linear aspH and quadratic Psph artificial 
viscosity terms (usually ~ 1 and sometimes, in some specific 
cases, < 1) are chosen = 1 and = 2, respectively. In the 
viscous models, the viscous force contribution is represented 
by the divergence of the symmetric viscous stress tensor 
in the Navier-Stokes equation. A symmetric combination 
of the symmetric shear tensor times the particle velocity 
has been added to the energy equation as a viscous heating 
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contribution. Tiie bulk physical viscosity contribution has 
not been considered for the sake of simplicity. 

Artificial and turbulent physical viscosities are indepen- 
dent from each other. The artificial viscosity terms should 
be smaller than the physically viscous ones, otherwise the 
physical viscosity role would be negligible. The relevance 
of viscous forces could be even compa rable to the gas 
pressure forces, esp ecially if ass = 1 fLa nzafame et ahl 
I2OO6I : lLanzafam3l2009iV An analytical formulation, describ- 
ing thsjiumerical artific ial viscosity coefficient, is reported 
in iMolteni et al.l l|l99ll ): vsph = Cgh, where Cs is the 
sound velocity. According to such a definition, its ratio with 
the Shakura-Sunyaev viscosity coefficient vss = otssCsH 
is: vsph/vss = h/{assH), for each SPH particle. For 
h/H — 5 • 10~^, where H is the scale-height of the disc, 
Vsph /fss ~ 5- 10~^/ass. This implies that the role of arti- 
ficial viscosity could be significant, compared to the physical 
viscosity role, if small ass and large h values are taken into 
accou nt. According to lMurravl (| 19961 ) and to lOkazaki et al] 
(|2002l ') the shear viscosity vsph ~ O.lasPHCsh with vss = 
O-ss-artifCsH . Accordiug to their results, the numerical ar- 
tificial viscosity coefficient vsph is even smaller if aspH ~ 1. 
In fact, the ratio vsph/vss — 0-lh/{assH). Hence, for 
H/hr^ 10-^20, vsph/vss ^ {5-10-^ ^10'^) /ass- This im- 
plies that, the role of artificial viscosity can be comparable to 
the role of a very low physical viscosity, because of the corre- 
lation between the SPH artificial viscosity parameter aspH 
and the Shakura-Sunyaev viscosity parameter ass-aruf is: 
ass-artif ~ O.laspnh/H without any bulk viscosity con- 
tribution and supposing gas incompressibility (V • = 0). 
Notice that, according to these correlations, the Shakura- 
Sunyaev parameter ass-artif (non zero only for approach- 
ing particles) is not the Shakura-Sunyaev viscous parameter 
ass for physically viscid gases, but the transformation of 
the artificial viscosity term into the Shakura-Sunyaev for- 
malism. Such results show that the gas compressibility has a 
relevant role since the physical viscosity mainly works when 
the density varies on a length-scale of the order of the ve- 
locity length-scale, not only as a bulk viscosity, but also as 
a shear viscosity. Moreover, notice that the assumption of 
an adaptive h SPH or a constant h SPH could also have 
a role both in artificial viscosity and in physical viscosity 
roles. These results show that the role of a fully viscous fiuid 
dynamics is still far from any conclusion and that physical 
assumptions as well as numerical hypotheses and boundary 
conditions are also determinant. 



3 VISCOUS FLUID DYNAMICS EQUATIONS 

As for viscous gas hydrodynamics, the relevant equations to 
our model are: 

-f pV • u = 

at 



continuity equation(l) 
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energy equation(3) 
perfect gas equation(4) 

kinematic equation(5) 



The most of the adopted symbols have the usual mean- 
ing: d/dt stands for the Lagrangian derivative, p is the gas 
density, e is the thermal energy per unit mass, '^grav is the 
effective gravitational potential generated by the two stars 
and uj is the angular velocity of the rotating reference frame, 
corresponding to the rotational period of the binary system. 
Self-gravitation has not been included, as it appears irrele- 
vant. The adiabatic index 7 has the meaning of a numerical 
parameter whose value lies in the range between 1 and 5/3, 
in principle, r is the viscous stress tensor, whose presence 
modifies the Euler equations for a non viscous fluid dynam- 
ics in the viscous Navier-Stokes equations. 



CLASSICAL SPH KERNEL AND PARTICLE 
SMOOTHING RESOLUTION LENGTH 



In its original formulation iMonaghanl 1 19851 . Il992l : 
iMonaghan fc Lattanziol 119851 ) Gaussian Kernels Wcij — 
W{rij,h) as: 



in ID 



(6) 



% 1 • 



, in 3D 



have been adopted in SPH, where Vij = \rij \ = \ri — Tj ] 
represents the module of the radial distance between par- 
ticles i_aBd_2;_Also, an example of "Super Gaussian Ker- 
nel" (|Monaghanlll992l ) has also been described. Even a fac- 
torization of Gaussian Kernels for each dimension h as also 
been adopted jShapiro et al.lll996l :[ Owen et al.lll998l ) in an 
ASPH formulation, adopting 3D ellipsoid Kernel geometry 
to achieve a higher accuracy, according to an anisotropic 
Vp-dependent spatial particle concentrations; or accord- 
ing to the mean particle spacing, as it varies in time , 
space, and direction around each particle (|Liu et al.ll2006l ). 
Kernels based on cubic splin e s since the end of the 80's 
(|Monaghan fc Lattanzicl 1 19851 : iMonaghanI Il992l ) have also 
widely been adopted. Typically, in 3D, such cubic spline 
Kernels W{rij, h) are in the form: 



+ if s; q^, ^ 1 







if 1 ^ qij ^ 2 
otherwise, 



(7) 



where qij — rij/h. 



5 GASPHER AN ALTERNATIVE WAY FOR 
KERNEL AND SMOOTHING LENGTH 

5.1 Lack of SPH and ASPH in FE conditions 

The hidden problem is whether ASPH, as previously formu- 
lated, are effective whatever is the compressibility regime 
considered, especially when FE conditions are adopted on 
the edges of the particle envelope. High compressibility gas 
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dynamics prevents us from distinguishing the truth regard- 
ing whether a SPH-like technique is able to perform a correct 
fluid dynamics, since accretion discs would be formed any- 
way even in physically inviscid conditions. In this case, the 
roles of the Kernel choice and of its resolving power are hid- 
den. Gas loss effects in low compressibility conditions natu- 
rally develop, especially at the disc's outer edge, because of 
the pushing action towards the outer space of particles just 
below the disc's surfaces and below the disc's outer edge, if 
the gravitational field is not able to keep gas particles in the 
gravitational potential well. In ASPH interpolation particle 
domains swell at both free edges (inner and outer). Nor- 
mally, in an accretion disc, the density is a decreasing func- 
tion of the radial distance from the central star. This implies 
that particle adaptive h should decrease towards the inner 
disc bulk, without any restriction imposed on the number 
of particle neighbours. Problems deriving from the inade- 
quacy of artificial viscosity role and the particle interpola- 
tion/interpenetration could be relevant. Even the choice of 
a threshold value for hmin as a lower limit would be arbi- 
trary and no differences would appear in results compared 
to classical SPH results, adopting the same h — hmin- If 
ASPH is adopted, even restricting the particle neighbours 
to a fixed number in its conservative form, the behaviour of 
h for each particle is contrary, swelling also within the disc 
bulk and producing enhanced gas loss effects at the disc's 
outer edge and on the disc surfaces, in spite of the viscos- 
ity eventually introduced, as well as a draining effect of the 
disc's inner edge toward the central compact star due to a 
stressed radial transport. 

Whenever and wherever spatial isotropy and homogene- 
ity hold, a modulation of spatial smoothing resolution length 
does not affect results, in principle, in so far as h is large 
enough to prevent particle interpenetration and neighbour 
particles are enough to allow good interpolations. However, 
the situation is rather different if spatial gradients exist. 

It is quite normal that a smaller threshold limit hmin 
is imposed on particle h because problems on the ineffec- 
tiveness of artificial viscosity in handling shock fronts would 
arise if /i — >■ 0, together with a too short time step com- 
puted according to the Friedrich-Courant-Lewy conditions. 
Artificial viscosity vanishes when the limit value of the par- 
ticle interpolation domain goes to zero, due to the fact that 
in its analytical expression it is linearly dependent on the 
smoothing length h. Its role, limited to a filing effect, should 
not be dominant compared to gas pressure terms. Therefore 
such condition is fully altered if, according to eqs. (18, 19), 
rjij ~ 1, see App. A. Other formulations of artificial viscosity , 
depending on particle mutual distance rij (jMonaehanll 19971 ) . 
do not modify the problem. Therefore, ASPH results would 
be deeply infiuenced by a dominant role of artificial terms 
if such a condition is mostly realized in sonic and subsonic 
regimes and/or in progre ssive turbulent rarefying regi mes . 
Some authors (|Morris fc Monaghanlll997l : [Owen et al.lllQQ^ ') 
handle artificial viscosity switching it off, especially in low 
density conditions, when particle h increases or in high tem- 
perature conditions when particle sound velocity is subsonic. 
As a result, the switching on/off of the artificial viscosity 
limits its role, but low density sonic and subsonic conditions 
stay still be critical. 

In both situations, the problem of a correct hydrody- 
namics involves not only the bulk of the gas structure in the 



computational domain, but mainly the physics of the FE 
of the computational domain. In particular the outer one 
for gas expansion problems and the inner one for collapse 
problems. 

As for physically viscous ASPH simulations, mass and 
angular momentum transport are deeply affected by the h 
particle smoothing resolution length. We expect a higher 
particle transport when particle h statistically increases and 
the opposite effect wh en h statistically decrea s es. In a low 
compr essibility regime. [Lanzafame et al.l l|2006l ) : lLanzafani3 
showed that physical turbulent viscosity hampers 
particle repulsion, due to pressure forces, contributing to 
accretion disc consistency and limiting particle loss at the 
disc's outer edge. However, if particle smoothing resolution 
length h increases in ASPH, and radial transport becomes 
unnaturally too much effective, the opposite effect arises so 
much that the inner edge of the disc could be indefinite. 



5.2 The Kernel of GASPHER: comparison to 
other Kernels 

In GASPHER modelling, a radial Gaussian-derived Kernel, 
related to the well-known "Error Function" with a constant 
smoothing length h equivalent to its HWHM is considered: 



WErF,ij 



1 

1 

27r3/2hr2. 



in ID 
in 2D 
in 3D. 



(8) 



In such a Kernel we stress that its interpolation ra- 
dial extension is unlimited, although its typical smoothing 
length h is spatially and permanently constant. In GAS- 
PHER, to collect an adequate particle neighbours number 
is not a problem because of the unlimited spatial extension 
of its Kernel. In the continuum limit, the three interpolation 
Kernels give the same interpolation integrals for ID flows, 
as well as the last two Kernels give the same interpolation 
integrals for 2D flows. 

The origin of this Kernel function relies in the well 
known "Error Function": 



ErF{x) 



e dt, 



whose "Complementary Error Function" is: 
2 



ErFC(x) = 1 - ErF{x) 

For X — 0, 
ErFC{0) = 1 - ErF{0) 



e dt. 



e"* dt = 1. 



(9) 



(10) 



(11) 



For a:: = 0, ErFCiQ) equals the zero order Gaussian 
integral: 



h = 



e dt=- — . 



(12) 



In performing 3D integral. 



WErF,ijd Ti. 
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1 2 2 

1 2 2 

2 f°° - 2 
= -r= / e qij=rij/h. (13) 

Hence, J WErF,ijd^rij = 1. 

Also J H''3s.ijrf''*rij — 1, as well as J Wcijd'^rij — 1, 
this last, considering the well known properties of Gaus- 
sian integrals: /„ = t"e~^* dt, and in particular I2 = 

J^t^e-'^dt^Tv^/^A. 

Fig. 1 displays, Wss.ij Wc.ij and Wetfaj as a function 
of qij = nj/h. AnqfjW3S,ij , AirQijWcij and 'iTvqfjWErF,ij, as 
well as Airqfj'VWss.ij , A-wqfj'VWcij and ■inqfj'VWErF^ij are 
significant for 3D integrations. Fig. 1 displays the much bet- 
ter GASPHER interpolation capabilities, with respect to the 
current SPH or ASPH techniques using other Kernels, not 
only because 3D interpolations are more weighted toward 
rij — > 0, but also because —VWErF,ij 00 as it should be, 
avoiding the well known "particle pairing instability" effect, 
affecting the other two behaviours (VWij displays a mini- 
mum for Tij/h « 1). In the conversion from mathematical 
integrals to computational summations in 3D, the role of 
A-Krfjd/rij is equivalent to . Thus, wherever ^ , 
and spatial gradients exist, the effectiveness of the adopted 
interpolation Kernel comes out. In the resolution of the Eu- 
ler or of the Navier Stokes equations, spatial derivatives have 
to be calculated. In the calculation of Vp/p in the momen- 
tum equation, two particles cannot coincide because the 
pressure force is physically infinite. Moreover, also for the 
V • u in the energy equation or in the continuity equation, 
this non physical case should be carefully avoided because 
no velocity divergence can exist if particle mutual separation 
is zero. Summing up, both indexes i 7^ j and rij > 0. In the 
unrealistic case of = 0, spatial derivatives to compute 
gradients or divergences can be bypassed because unphysi- 
cal. In the case of a very short particle mutual separation, 
natural computational difficulties can arise only for a very 
small particle separation. This is unavoidable when a very 
high compression characterize the fluid, because calculated 
pressure and individual pressure forces are always naturally 
very high. However, in particular for accretion or collapse 
processes, the particle merging in a new particle, created at 
the centre of mass, conserving mass, energy and momentum 
could be the best solution. This is a useful physical expedi- 
ent, also used in ASPH, whenever a strong gas compression 
occurs. It avoids a too short explicit time step calculation, 
according to the well known Fricdrich-Courant-Lewy. In par- 
ticular, for ASPH technique only, it also avoids any artificial 
viscosity inadequacy in handling shocks. Such an allowed ex- 
pedient could be correctly also used in GASPHER in such 
conditions. 

We pay attention that in 3D interpolations, it is not the 
role of the Kernel W that it is important. Instead, it i.s the 
Wr^ that is to be taken into account, as Fig. 1 clearly dis- 
plays. Hence, on ith particle, when nj — >■ 0, Wr^ converges 
toward a finite value. If this is the explanation regarding 
the continuum limit, in the spatial discretization, this role 
is carried out by the particle density pj which, in the SPH 
formulation Ai = '^^AjW^j/nj = '^i'mjAjWij/pj, divides 
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Figure 1. Radial plots of SPH Kernels W^Sjij (^q- 7), Wa,ij 
(cq. 6), as well as of GASPHER WErF,ij (eq. 8). qij = Tij/h. 
'iTrq?.W3s,ij, "iTrqfjWcij, '^■Kq'ijWErF ,ij, as well as the radial 
derivatives times 47rq|^ both useful for 3D calculations are also 
reported. 

Wij. Only when — >• GASPHER formulation becomes: 
Ai = ^-AjWij/uj AnAjWijVijh^. On the other hand, 
whenever nj — > the concept of dimension is meaningless. 
This implies that, if rtj — ^ 0, and especially if rij = 0, the 
ID formulation of Kernel can also be taken into account to 
simplify computational complications in some selected cases, 
whenever the 2D or the 3D fluid kinematics flows along one 
selected direction. 

5.3 Advantages of GASPHER 

This Kernel choice resolves the problem of neighbours in- 
adequacy, as well as the problem of the SPH and ASPH 
"particle pairing instability" for rij/h < 1 due to the fact 
that when nj 0, — Vp does not become infinite. 

For practical reasons in computational resources, even 
a limitation to several h of the order of Ih with ? ~ 4 10 
could be considered with very small modifications in results, 
keeping constant the resolving power h of all particles. In 
fact, theoretically considering a homogeneous and isotropic 
3D particle distribution, if Annh^ /3 (where n is the parti- 
cle concentration) is the number of neighbours closer than h 
for each ith particle, it increases up to 4-Kn{lh)^ /3 i.e. up to 
~ 64 ^ 10^ times. Alternatively, neighbours can be limited 
to a selected number (40 in our 3D models). In both cases, 
a very small modifications in results, neglecting further in- 
terpolating particles, is made because the most important 
neighbours in the interpolation are the closest ones. In this 
case, if neighbours are a large number, due to a very high 
particle concentration, it is easily possible to merge more 
particles in a single new particle, created at the centre of 
mass, conserving mass, energy and momentum. So doing, 
the ASPH's risk to decrease the spatial smoothing resolu- 
tion length to values involving an ineffective artificial vis- 
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cosity behaviour, as well as the danger to get a too small 
computed time step in the Courant-Friedrich-Lewy condi- 
tion when /i — > 0, are avoided. A fixed number of neigh- 
bours can be a serious risk by limiting the interaction to 
30 -T- 40 neighbours only the inner "flat part" of the Kernel 
contributes to SPH sums. A large contribution from other 
particles outside this "flat part" would be wrongly neglected. 
In GASPHER this does not occur because the Kernel slope is 
not "flat" for rij — >■ 0, instead VWErF.ij — > — oo. Of course, 
also ASPH techniques try to avoid the unpleasant "particle 
pairing instability". However, the particle resolving power 
h cannot decrease too much in regions of very high parti- 
cle concentrations otherwise the artificial dissipation due to 
the artificial viscosity does not work well. Moreover, at the 
same time, the time step explicitly computed according to 
the Friedrich-Courant-Lewy condition becomes too short if 

The possibility of adoptin g a numerical SPH code, in- 
cluding the physical viscosity (|Flebbe et al.l 



sidering lLanzafame et al] (|2006l ^ jLanzafarag 

results. 



(2008a!|a, 120091) 

makes us able to answer the problem whether 
ASPH's and/or GASPHER methods are reliable in improv- 
ing fluid dynamics compared to the original SPH, where 
the s moothing length h is constant. Although some au- 



thors l|Fulbright et al.ll 19951 : [Shapiro et al.lll99d : lQwen et all 
1 19981 ) adopted Gaussian Kernels, their methods belong to 
the ASPH numerical schemes where a spatially and tem- 
porarily variable smoothing length h is adopted. Although 
many efforts try to conciliate a reliable adaptive interpola- 
tion technique with computational resources, ASPH meth- 
ods are unsatisfactory in describing a correct gas dynamics 
because hidden numerical errors exist inside an adaptive in- 
terpolation, better revealed in a viscous transport process in- 
side a definite potential well. All ASPH's difficulties in han- 
dling the artificial viscosity dominant role in subsonic and/or 
expanding regimes, as discussed before, are prevented in 
GASPHER by the fact that the particle resolving power 
h is constant and equal to HWHM of spatially unlimited 
Gaussian Kernels. GASPHER technique limits the problem 
of particle disorder in compu ting p article V -v, as dis cussed 
in llmaeda fc Inutsukal (|2002l l and in lMonaehanl ()2006l ) as far 
as shear fiows are concerned because, even considering dis- 
ordered flows, particle disorder is tamed by GASPHER ex- 
tended interacting particle domains. In fact, the longer the 
particle interpolation range, the better the computational 
result, without any modification of particle resolving power 
h. 

Finally, an adequate fixed smoothing resolution length 
h allows us to resolve gas turbulence within the confined in- 
tegration domain even in low compressibility regimes. In non 
viscous conditions the local Reynolds n umber Re = Lv/v ^ 
IQ^v/cs, considering L = 0.5, v « Csh ijMolteni et al.lll99ll l 
and, more str e ssing, Re ^ IQ^v/cs considering v k, Q.lCsh 
l|Murravlll99i lOkazaki et alll2002l ). because of aspH = 1. 
Being the whole disc structure typically supersonic, even for 
7 — 5/3, Re > 10^. Hence, a moderate turbulence is ef- 
fective in non viscous 7 conditions, where gas collisions are 
relevant. Instead in viscous conditions, in the Shakura and 
Sunyaev formulations, no turbulence is recorded. 

If an adaptive method is adopted in low compressibil- 
ity conditions, the increasing of the smoothing resolution 
length h, up to an order of magnitude, prevents any tur- 



bulence resolution in an accretion disc, even for supersonic 
regimes. The evaluation of the minimum linear dimensions 
of the integration domain, able to solve turbulence adopting 
an ass parameter of the order of 0.1 — 0.5, gives a value of 
the order D ^ 10^^ ~ 10^^ in order to get a Reynolds num- 
ber Re ^ 10^, the smaller value is for supersonic regimes. 
The integration domain (the length of the primary's poten- 
tial well) of the order of 0.5. Therefore, how to handle an 
adaptive SPH with the problem of solving the turbulence is 
a real difficulty, and the adopted fixed h is correct in order 
to solve this problem. Larger (and adaptive) h values are 
in open conflict with D ^ 10~^ -t- 10~^ in order to solve 
turbulence. 

These conclusions on turbulence in accretion discs hav- 
ing free edge boundaries are not those concerning the con- 
cept of turbulence wherever fixed static boundaried are con- 
sidered. Whenever particles move within a confined box, 
both SPH and ASPH results are traditionally correct in so 
far as Vij/h is not too small. In this case the problem regards 
the particle chaotic collisions in a close environment where 
the particle mean free path is less than two or three times 
the particle smoothing resolution length. 



6 GASPHER DISC SIMULATIONS IN CB: 
RESULTS AND DISCUSSION 

Looking at our SPH r esults in a physically viscous low 
compr es sibilit y regime (iLanzafame et ahl l2006l : iLanzafamd 
l2008al lbl [2009l ) as a reference, where the particle smoothing 
length is constant and a typical cubic spline function as a 
smoothing function have been assumed, we systematically 
perform a series of GASPHER simulations with the aim of 
getting a physically viscous well-bound accretion disc in a 
close binary. We show that such transport phenomenology, 
in a low compressibility regime, is significant in deciding the 
reliability of the adopted Kernel formulation for SPH fluid 
dynamics simulations, especially whenever free edge bound- 
ary conditions must be taken into account. 



6.1 Parameters and boundary conditions 

The characteristics of the binary system are determined by 
the masses of the two companion stars and their separa- 
tion. We chose to model a system in which the mass Mi of 
the primary compact star and the mass M2 of the secondary 
normal star are equal to IMq and their mutual separation is 
di2 = lO*' Km. The primary's potential well is totally empty 
at the beginning of each simulation at time T = 0. The injec- 



tion gas velocity at LI is flxed to Wi, 



130 Km while 



the injection gas temperature at LI is fixed to To = 10'' K, 
taking into account, as a first approximation, the radiative 
heating of the secondary surface due to lightening of the 
disc. Gas compressibility is flxed by the adiabatic index 
7 = 5/3. Supersonic kinematic c onditions a t LI a re dis- 
cussed in ILanzafame et al.l l|2006l ): lLanzafani3 l|2009l ). espe- 
cially when active phases of CB's are considered. However, 
results of this paper are to be considered as a useful test to 
check whether disc structures (viscous and non) show the 
expected behaviour. The reference frame is that centred on 
the primary compact star and corotating, whose rotational 
period, normalized to 2-k, coincide with the orbital period 
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of the binary system. This explain why in the momentum 
equation (eq. 2), we also include the Coriolis and the cen- 
trifugal accelerations. 

In our models the unknowns are: pressure, density, tem- 
perature, velocity, therefore we solve the continuity, mo- 
mentum, energy, and state (perfect gas) equations. In or- 
der to make our equations dimensionless, we adopt the fol- 
lowing normalization factors: = Mi -I- M2 for masses, 
di2 = 10" cm for lengths, Va = (G(Mi M2)/di2)^/^ 
for speeds, so that the orbital period is normalized to 2-k, 
po = 10~^ g cm~^ for the density, po = povl dyn cm~^ 
for pressure, vl for thermal energy per unit mass and To = 
(7 — l)v1 rup Kg^ for temperature, where rup is the proton 
mass and Kb is the Boltzman constant. The adopted Kernel 
resolving power in the GASPHER modelling is ft = 5 ■ 10~^. 
The geometric domain, including moving disc particles, is a 
sphere of radius 0.6, centred on the primary. The rotating 
reference frame is centred on the compact primary and its 
rotational period equals the orbital one. We simulated the 
physical conditions at the inner and at the outer edges as 
follows: 

a) inner edge: 

the free inflow condition is realized by eliminating particles 
flowing inside the sphere of radius 2 • 10~^, centred on the 
primary. Although disc structure and dynamics are altered 
near the inner edge, these alterations are relatively small 
because they are counterbalanced by a high particle con- 
centration close to the inner edge in supersonic injection 
models. 

b) outer edge: 

the injection of "new" particles from LI towards the inte- 
rior of the primary Roche Lobe is simulated by generat- 
ing them in fixed points, called "injectors", symmetrically 
placed within an angle having LI as a vertex and an aper- 
ture of ~ 57°. Normally, as a dopted since our fir st paper 
on SPH accretion disc in CB (jMolteni et al.lll99ll ). the ra- 
dial elongation of the whole ensemble of injectors is ~ 10ft. 
The initial injection particle velocity is radial with respect 
to LI. In order to simulate a constant and smooth gas injec- 
tion, a "new" particle is generated in the injectors when- 
ever "old" particles leave an injector free, inside a small 
sphere with radius hmin, centred on the injector itself. Par- 
ticle masses are determined by the assumed local density at 
the inner Lagrangian point LI: pLi = W~^g cm~^ (as typ- 
ical stellar atmospheric value for the secondary star), equal 
to m = pLi{hdi2f/{Mi + M2). 

The formulation adopted for the 3D SPH v i scous accre- 
tion disc models is the well - know n ggg lShakural (|l972l . [l97^ ) 
and IShakura fc SunvaevI 1 19731 ) parametrization: uss = 
assCsH, where Cs is the sound velocity, ^ ass 5; 1 and 
H — T^yCsliMxjr^yY^'^ is a dimensionless estimate of the 
Standard disc thickness, where Vxy — (Xl -I- Y'^'')^^'^ is the 
cylindrical radial coordinate of the ith particle. In this paper 
we adopt ass = 1 to point out evident differences in disc 
structure and dynamics between our disc models. 



6.2 General results 

We carried out our low compressibility (7 = 5/3) simula- 
tions until we achieved fully stationary configurations. This 
means that particles injected into the primary potential well 
(which is not deep, according to the primary small mass) are 
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Figure 2. XY plots and txy^ plots for both the inviscid (ass = 
0) and the viscous (oss = 1) disc model. The final time T and 
the total particle number A^, as well as the injection velocity from 
the inner Lagrangian point LI and ft, are also reported. 



statistically balanced by particles accreted onto the primary 
and by particles ejected from the outer disc edge. 

The orders of magnitude of the mass transfer injection 



rate from LI: Minj, the accretion rate Mace 
rate Meje are « lO", « 5.5-10^'' and « 4.5 
viscous model and ~ 8.0 ■ 10^ 
the viscous model, respectively. poh?'d\2Vo, is the conversion 



and the ejection 
10^®, for the non 
7.9 • lO^'' and -10^^ for 



fa ctor from particle/time t o ff s 
in lLanzafame et ahl \ 



Such values (al so ado pted 

20061 ) ; iLanzafamd (|2008al lbl, bOOgl ) ) are 

representative of active phases of CB whenever either the 
restricted problem of three bodies in terms of the Jacobi 
constant or the Bern oulli's theorem are t aken into account 
during such phases (|Lubow fc Shul Il975l ) , considering the 
conservation of the flux momentum in the crossing of LI 
from the two Roche lobes. 

Fig. 2 displays XY plots of both the physically invis- 
cid and viscous disc models (ass = and ass = 1). N 
represents the total number of particles in each model. In 
classical SPH, no well-bound structures with a definite disc's 
outer edge come out (Molteni et al. 1991; Lanzafame et al. 
1992, for 7 = 1.1 and 7 = 1.2). The inviscid GASPHER disc 
model shows a higher particle concentration at the disc's in- 
ner edge, close to the primary star. Instead, a well-defined 
structure comes out in the viscous disc model in stationary 
conditions. 

Fig. 2 also displays the rxvZ plots obtained by folding 
all disc bulk particles onto a plane containing the Z axis and 
being perpendicular to the XY orbital plane. An evident lati- 
tudinal spread appears for the inviscid model. Computed lat- 
itudinal angular spread is « 60° fo r inviscid disc model. This 
res ults compare to that ob tained in iMolteni et al] l|l99ll ) and 
in iLanzafame et aLlll992l). as far as the non vi s cous model, 
and to that obtained in Lanzafame et al.1 ()2006h : iLanzafarng 
(|2009l ) as far as the viscous model, are concerned. As for the 
viscous model, the latitudinal spread is ~ 24°. 
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Fig. 2 clearly displays the coming out of spiral patterns 
in the XY plot of the viscous m odel. These particular st ruc- 
tures did not come out in SPH iLanzafame et al.l (|2006') re- 
sults, where both the same supersonic injection conditions 
from LI, and the same stellar masses, as well as the same 
low gas compressibility were adopted in vis cid ass = 1 
cond i tions. However, an exhaustive lite r ature ( Sawada et al. 
1987'; 'S pruit et al * 1987; Kaisig 1989'; 'Sawada fc Matsudal 
1992: Savoniie et al. 1994: Lanzafarne et al,.. 2000 , .2001.) ex- 



ists, showing which conditions favour the development of 
such structures (e.g. tidal torques , external and/or outer 
edge perturbations). In particular, ILanzafame et all ()2000l . 
I2OOII ) showed that high angular momentum injection con- 
dition from LI produces these patterns. This beyond doubt 
shows the better effectiveness of GASPHER Kernel choice 
(33) compared to the common cubic spline SPH Kernel an- 
alytical formulation. 

The c ompar ison with ILanzafame et all l)2006l l: 
iLanzafamd (Hooi) results ensures us that GASPHER 
technique performs not only correct calculations, but also 
that particles at disc's outer edge are not isolated. Moreover, 
the full radial transport cannot be affected by any "particle 
pairing instability" because the Kernel formulation (33) 
prevents such unpleasant inconsistency in the disc bulk. 

Low compressibility gas loss effects affect the non vis- 
cous GASPHER disc su rfaces and out e r edg e. The same 
re sult were obtained in [ Molteni et al.l l|l99lf ). as well as 
in ILanzafame et all (|l992l ). working in SPH and adopting 
two different spatial smoothing resolution lengths and sonic 
injection transfer conditions from LI. Supersonic injection 
conditions fro m LI are now taken into account, as described 
m .Lanzafamel (|2009l ) for active phases of CB. Therefore, non 
viscous gas loss effects from disc's outer edge and surfaces 
are a fortiori correctly expected, taking into account of the 
higher injection mechanical energy from LI. At the same 
time, even though the low density non viscous disc struc- 
ture is statistically rarefied, no neighbour inadequacy affect 
GASPHER interpolation. 

The low total number of particles within the pri- 
mary's potential well (~ 3000) in non viscous conditions 
is due to the absence of any physical viscosity able to 
keep bound particles against pressure forces responsible of 
particle removal from the disc outer free edge for 7 — 
5/3 whenever low mass CB's are considered . This result 
is well known |Molteni et al.l Il99ll : ILanz afame et al.l 1 19931 : 
ILanzafamd l2008bl . l2009l ~ This particular is not trivial be- 
cause the bound of the edge of the computational domain 
prevents any particle removal allowing to get the wished par- 
ticle concentration in spite of the effective particle repulsion 
for high 7 values. 

6.3 Accuracy 

SPH free surface flows were developed bv lMonaghanI ()l994h . 
with the aim of solving the unpleasant problem of FE layer in 
SPH techniques, mainly to simulate breaking waves, but at 
relatively low resolution. A reduction in noise, with smooth- 
fre e surfaces and r egula r part icle distribut i on, wa s obtained 
by iBonet fc Lo5 (|l999l ) and iBonet et all (|2004 ). develop- 
ing SPH models where first order completeness was en- 
forced, that is that first order polynomials are exactly re- 
produced. Error estimates in a SPH interpolant are eval- 
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Figure 3. Logarithmic plots of specific angular momentum ra- 
dial distributions and of temperature for both the inviscid models 
{ctss = 0)1 S'IkI the viscous models (ctss = !)• The final time T 
and the total particle number N, as well as the injection velocity 
from the inner Lagrangian point LI and h, are also reported. 



uated in iMonaghanI () 19851 . \l99± . However in this paper, 
the lack of completeness of SPH interpolants is not taken 
into account. A formulation for the total error, determin- 
ing how simulation parameters should be chosen and taking 
into account of t he order of complet eness is still not written 
in the literature. iBonet et al.l l|2004 ) adopted modified Ker- 
nel gradients into the classical SPH equations. However, the 
hidden problem with this approach is that modified Kernels 
no longer have the property that spatial gradients with re- 
spect to their two position arguments are exactly opposite 
between two contact par ticles. This Kernel pr operty is es- 
sential in SPH equations. IVaughan et al.1 (I2OO8I ) showed " an 
expression for the error in an SPH estimate, accounting for 
completeness, an expression that applies to SPH generally" , 
paying attention to the conservation principles. They found 
that a common method, enforcing completeness, violates the 
conservation principle of Kernel spatial gradients must be 
opposite between two contact particles. They also showed 
some examples of discretization errors: numerical boundary 
layer errors. Errors for a SPH summation interpolant are 
functions o f both particle d i stribution and partic le smooth- 
ing length (|Monaghanlll985l : rVa"ughan et aLlboOSl ). In an ex- 
act formulation, such errors are described by both volume 
and surface integrals of both neighbour particle distribution 
and their smoothing resolution length h. Therefore, in FE 
layer conditions, not only relevant errors in interpolations, 
but also unnatural pressure gradients in FE conditions at 
the edge of the computational domain occur in ASPH. 

A reasonable SPH accuracy is related to the number of 
space neighbours of each SPH particle. As a free Lagrangian 
numerical method, classical SPH methods, are free from er- 
rors as far as momentum and angular momentum are con- 
cerned. Instead, errors can occur as to energy as for ASPH 
variants. In particular, it is remarkable the evaluation of the 
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energy error propagator (EEPR) {AE / dt) / E , computed for 
each particle, to have the correct idea of temporal propa- 
gation of energy errors. If SPH-like methods involve a sys- 
tematic error in energy [A_B/_E]|err of a few percent, this 
error progressively increases in time as J [{dE / dt) / E]\errdt. 
This means that, if a long time is necessary to achieve a fully 
stationary configuration, errors in energy conservation could 
be significant. The evaluation of the GASPHER EEPR for 
inviscid disc model is « 3.02-10"'^. Instead, the GASPHER 
EEPR for viscous disc model is ^ 1.85 • 10"^. Being errors 
in energy of this order of magnitude we do not usually al- 
low to distinguish if numerical simulations correspond to a 
fluid physical behaviour. Thus, the numerical error in energy 
on particle over-expansion/over-compression is not domi- 
nant step by step. Unfortunately, it accumulates in time. 
This implies that numerical simulations, limited only to ex- 
plosive or collapse short time tests, would not be reliable 
in testing ASPH codes. Therefore, once more, this conclu- 
sion strengthens numerical tests and simulations based on a 
transport mechanism. 
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6.4 The role of physical viscosity 

Physical viscosity naturally works where the particle mutual 
velocities (and separations) change in time, namely when a 
mutual acceleration exists, contrasting gas dynamics (rar- 
efaction or compression) and converting kinetic energy in 
thermal energy. Such a mechanism clearly supports the de- 
velopment of well-bound accretion discs inside the primary 
potential well, in spite of the low compressibility, at least for 
ass = 1 both in classical SPH and in GASPHER approach. 

We want to point out that adopting ass = 1 does help 
emphasize differences in disc structure and dynamics com- 
pared to the physically inviscid model. However, values of 
ass smaller than the unity may be m ore realistic accord- 
ing to some thin disc analytical models (jPringle et al.|[l98^ : 
lLa8otall200ll l. We recall that our physical viscosity is only a 
shear viscosity. For the sake of simplicity, no bulk viscosity 
has been considered, as explicitly mentioned in the paper. 
In fact, a value 035 = 1 for the bulk viscosity should be 
too high. Fig. 3 displays, in a logarithmic scale, the angu- 
lar momentum and temperature radial distribution, for all 
models. Such radial distributions for the GASPHER vis- 
cous model are very close to that of the Standard model 
r^il cx r^^^ and T oc r""^''*. This can be explained consid- 
ering that, in stationary conditions, an accretion disc redis- 
tributes the angular momentum injected at the outer edge 
into the disc bulk, according to outer edge boundary c ondi- 
tions only, as already s hown in [Belvedere et al] ()l993l ) and 
iLanzafame et all l| 19931 ) . Physical viscosity plays a role in re- 
gions where particle velocity gradients are significant. This 
means that physical viscosity plays a relevant role mainly in 
the radial transport, while it has scarce influence on the tan- 
gential dynamics. A strong difference appears when looking 
at the temperature radial distribution. In fact, the heating 
effect of the physical viscosity is particularly evident in the 
disc's inner zones. We recall that the disc itself is in an 
equilibrium stationary state where the heated particles are 
directly accreted towards the primary. This as far as parti- 
cle advection is concerned. As for conduction, although it is 
much less important, notice that the temperature decreases 
towards the exterior, thus dispersing heat outside. However, 



Figure 4. XY plots and rxY^ plots for the inviscid and the 
viscous [ass = and ass = 1) disc models for a different particle 
smoothing resolution lengths h = A ■ 10~^. The final time T and 
the total particle number N, as well as the injection velocity from 
the inner Lagrangian point LI, are also reported. 

discs could also radiate energy. In disc models without ex- 
plicit inclusion of radiative terms in the energy equation 
(almost all models, since, this inclusion complicates things 
considerably), the effect of radiative cooling is better simu- 
lated with 7's less than 5/3. 



7 INFLUENCE OF THE GASPHER 

SMOOTHING RESOLUTION LENGTH ON 
DISC STRUCTURE: DOES THE SPATIAL 
RESOLVING POWER AFFECT GENERAL 
RESULT? 

To study how EE fluid dynamics is affected by the initial 
smoothing resolution length choice, we performed two more 
simulations (non viscous and viscous) adopting a smaller 
smoothing resolution length: h — i ■ 10~^, improving spa- 
tial resolution since particle injection. Taking into account 
of injection condition from our previous simulations for 
hi^Li = 5 • 10~^, particle masses are scaled, conserving 
the same mass density from LI, according to the ratio of 
particle volumes: rrii^h = rrti^h^o.oos (/ii,Li/5- 10"^)^ Thus, 
the mass transfer rate from LI is self-consistent and auto- 
matically comparable to that relative to simulations with 
hi^Li ~ 5 • 10~^, without any variation of the injection ve- 
locity Vinj ~ 130 Km . To do this, it is necessary to 
recalculate newly the total number of injectors, by adopting 
the simple scale law: Ni„j = A''i„j,h=o.oo5(5 • 10~^//ii,Li)^. 
Hence, according to these simple scaling laws, we keep in- 
jection conditions comparable both for the initial density 
and for the mass transfer rate at LI. 

Results of such further simulations are displayed in 
Figg. 4 and 5, where plots of such 7 = 5/3 3D GASPHER 
simulations are displayed, free of any difficulty on the suf- 
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Figure 5. Logarithmic plots of the radial distribution of the spe- 
cific angular momentum and of temperature for both viscous and 
non viscous accretion discs. «ss, the particle smoothing reso- 
lution lengths h, as well as the final time T, the total particle 
number N, and the injection velocity from the inner Lagrangian 
point LI, are also reported. 



ficient number of neighbour particles, as explained before. 
The total number of disc particles shows a monotonia in- 
crease by decreasing h. Thus, both disc density is compara- 
ble in both disc models as well as the mass of the simulated 
discs being Nh^, as well as Nrrii ~ constant, within statis- 
tical fluctuations. 

In the non viscous regime the larger number of disc 
particles are still affected by a gas chaotic coUisional com- 
ponent on top of the spiral disc's kinematics. Moreover, the 
Reynolds number increases because of the reduction of the 
particle smoothing resolution length, from h = 5 ■ 10~^ to 
/i = 4 • 10"^. Instead, whichever is the GASPHER adopted 
particle smoothing resolution length h, in a viscous regime, 
both the radial transport of mass and angular momentum, 
as well as the radial temperature profile, are not sensitive to 
any adopted particle resolving power as Figg. 2 to 5 clearly 
display. Their radial behaviour is strictly comparable to that 
of the typical standard disc, whose specific angular momen- 
tum r^n oc r^/^ and whose mean temperature T oc r"^''^. 
Hence, this result is a further confirming check that GAS- 
PHER result, in their general aspect, as far as the radial 
transport and thermal properties, are not strongly depen- 
dent on the assumed spatial resolution. Moreover, the local 
physical properties are clearly comparable with each other, 
being the particle spatial resolution in the graphs different, 
but not their physical values, that is denser (more rarefied), 
lighter (heavier) particles to get the same density, as an ex- 
ample. 



8 CONCLUDING REMARKS 

From the astrophysical point of view, our results show that 
in GASPHER modelling, where particle interpolation ra- 



dial extension is conceptually unlimited - although parti- 
cle smoothing length h is spatially and permanently con- 
stant - solve the problem of neighbours inadequacy. More- 
over, physical viscosity supports the development of a well- 
bound accretion disc in the primary potential well, even in 
the case of a l ow compressibility gas dynamics. Such resul ts, 
also shown in lLanzafame et al.l (|2006l ): lLanzafamel (|2008al lbl. 
120091), mean, once more, that the initial angular momentum 
injection conditions at the disc's outer edge are responsible 
for the disc tangential dynamics, while viscosity is mainly 
responsible for the thermodynamical disc properties, even 
for low compressibility disc models (7 > 1.1, here 7 — 5/3) 
when gas loss effects are physically expected according to 
the low compressibility gas dynamics and to the low stellar 
mass of the central accretor. Moreover, in GASPHER vis- 
cous fluid dynamics, further details of the flow are revealed 
(e. g. the coming out of spiral patterns in disc structures). 

From the numerical point of view, reliable results are 
reproduced in a GASPHER, despite EE conditions are 
adopted. Without considering the injected particle stream, 
such simulations could also be considered as accretion and 
transport general tests within a gravitational potential well. 
Typical tests as far as non viscous ID shock tube show that 
GASPHER technique produce results in a very good com- 
parison with analytical ones, having the advantage to solve 
the EE difficulties without any "particle pairing instability". 
Simulation, carried out in low compressibility and in high 
viscosity conditions, to stress out results, is significant to un- 
derstand the quality of numerical code. The transformation 
of SPH codes in a GASPHER code, without further numeri- 
cal efforts, seems likely to be an interesting future challenge. 
As far as the computational cpu time is concerned, there is 
not conceptually any disadvantage in such transformation, if 
particle neighbours are fixed (e.g. 30 or 50) for each particle, 
by the introduction of a boundaries counter/limiter because 
the number of particle neighbours rules the computational 
cpu time. 

The necessity to perform better SPH numerical inter- 
polations on contact surfaces, or at EE layers, recently in- 
spired authors to develop SPH-derived techniques to achieve 
a higher accuracy. An SPH dynamic refinemen t has re- 
cently been developed bv lFeldman fc Boned l|200'ii) to calcu- 
late boundary contact forces in fluid flow problems through 
boundary particle splitting. Such a technique could also be 
very interesting and competitive in solving EE problems. 
However, this is beyond the scope of this paper. 

We conclude that although high compressibility invis- 
cid results among different schemes could compare with each 
other especially if constraints are imposed on boundaries 
of the computational domain, differences arise either if EE 
and/or if viscous flows are involved. In such conditions, GA- 
SPHER technique shows a regular behaviour and better con- 
serve the total energy, as well as reduces the influence of 
the artificial viscosity for non viscous ideal shear flows free 
of any gas compression (see Appendix). Computational cpu 
time is mainly governed by the number of neighbour parti- 
cles for each particle. Therefore, no disadvantages arise, in 
principle, in adopting a GASPHER code with respect to an 
ASPH code if the neighbour particle statistical number is 
the same. 
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APPENDIX A: SPH FORMULATION OF BOTH 
PHYSICALLY INVISCID AND VISCOUS 
PERFECT GAS HYDRODYNAMICS 

Al SPH and ASPH (in adaptive smoothing 
length h) techniques 

The SPH method is a Lagrangian scheme that discretizes 
the fluid into moving interacting and interpolating domains 
called "particles". All particles move according to pressure 
and body forces. The method makes use of a Kernel W use- 
ful to interpolate a physical quantity A{r) related to a gas 
particle at position r according to: 



A{v) 



A{r')W(r,r' ,h)dr' 



(Al) 



W{r,r',h), the interpolation Kernel, is a continuous 
function - or two connecting continuous functions whose 
derivatives are continuous even at the connecting point - 
defined in the spatial range 2h, whose limit for /i — > is the 
Dirac delta distribution function. All physical quantities are 
described as extensive properties smoothly distributed in 
space and computed by interpolation at r. In SPH terms we 
write: 



A, 



spAj_ 



Wiri,rj,h) 



(A2) 



where the sum is extended to all particles included 
within the domain D, rij — pj/rrij is the number density 
relative to the jth particle. W{ri,rj,h) ^ 1 is the adopted 
interpolation Kernel whose value is determined by the rela- 
tive distance between particles i and j. 

In SPH conversion of mathematical equations (eq. 1 to 
eq. 4) there are two principles embedded. Each SPH parti- 
cle is an extended, spherically symmetric domain where any 
physical quantity / has a density profile fW{ri,rj,h) = 
fW{\ri — rj\, h) = /W(|r,;j|, h). Besides, the fluid quantity 
/ at the position of each SPH particle could be interpreted 
by filtering the particle data for f{r) with a single window- 
ing function whose width is h. So doing, fiuid data are con- 
sidered isotropically smoothed all around each particle along 
a length scale h. Therefore, according to such two concepts, 
the SPH value of the physical quantity / is both the overlap- 
ping of extended profiles of all particles and the overlapping 
of the closest smooth density profiles of /. This means that 
the compactness of the Kernel shape gives the principal con- 
tribution to the interpolation summation to each particle by 
itself and by its closest neighbours. In both approaches the 
mass is globally conserved because the total particle number 
is conserved. 

In SPH formalism, equations (2) and (3) take the form: 



dvi 
~dt 



E 



Pi f 



Pi 



+ 



(A3) 



E' 



+ 



E^ 



Pt 



Vi (Tj 

+ Vv3 — 



V,W^J (A4) 



where 



-2U3 X Vi + OJ X {UJ X Vi) ~ V<&9 



Vi — Vj , nij is the mass of jth particle and p* = Pi+ artificial 
pressure term. Ei = (e^ -I- ^vf). The viscous stress tensor r^/? 
includes the positive first and second viscosity coefficients rj„ 
and ^„ which are velocity independent and describe shear 
and tangential viscosity stresses (r/v), and compressibility 
stresses (Cd): 



Tafl = rfvCJaf} + (vV ' V 

where the shear 



— _|_ 

dXi3 dXa 



2 



(A5) 



(A6) 



In these equations a and /3 are spatial indexes while 
tensors are written in bold characters. For the sake of sim- 
plicity we assume ("i, — 0, however our code allows us also 
different choices. Defining 



Viafi 



N 

\ ^ rriiV 



1 ^jia 

dWij 



Pi dx0 



(A7) 



as the SPH formulation oidva/dxp, the SPH equivalent 
of the shear is: 



'iaP 



(A8) 



A full justification of this SPH formalism can be found 
in lFlebbe etal] l|l994al lb[). 

In this scheme the continuity equation takes the form: 



JV 

dpi _ \ ^ 
It ^ 2^ 



mjVij ■ ViWij 



or, as we adopt, it can be written as: 

JV 

^^mjWij 



(A9) 



(AlO) 



which identifies the natural space interpolation of par- 
ticle densities according to equation (9). 

The pressure term also includes t h e arti fi cial v iscos- 
ity contribution given by MonaghanI l| 19851 . 1 19921 ') and 
^Monagha n fc Lattanzid (| 19851 ). with an appropriate thermal 
diffusion term which reduces shock fluctuations. It is given 
by: 



llij = asPHfJ-ij + PsPHPij, 

where 



pij 



if Vi 



< 



(All) 



(A12) 







otherwise 



with Csi being the sound speed of the ith particle, 
^ /i^, aspH ~ 1 and /3s ph ~ 2. These aspH and 
PsPH parameters of the order of the unity are usually 
adopted to damp oscillations past high Mach number shock 
fronts developed by non-linear instabilities ()Boris fc Bookl 
197 3). These oisph an d Psph values were also adopted by 
iLattanzio et al.l l| 19851 ). Smaller aspH and Psph values, as 
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adopted bv lMeglicki et alj l|l993h . would develop more tur- 
bulence in the disc and possibly only one shock front at the 
impact zone between the infalling particle stream and the 
returning particle stream at the disc's outer edge. In the 
physically inviscid SPH gas dynamics, angular momentum 
transport is mainly due to the artificial viscosity included in 
the pressure terms as: 



2 ' 2 



+ 



Pi 



(A13) 



where p is the intrinsic gas pressure. 

The advantage of an ASPH is to perform better par- 
ticle interpolations ensuring a large enough nu mber of in- 
terpo lating particle neighbours. Several authors (|Benz et al.l 
I1990I ) have more recently adopted a criterion where the 
number of SPH particle neighbours for each time-step cal- 
culation is a fixed number, generally of the order of 30 
50, decoupling the h resolving power calculati on by any 
physical qua ntity. In stead, in previous p apers (Monaghan' 
1992; Fulbri ght et aLllTg gS: Sha piro et all [1996: Owen et al. 



19981 : iLiu et al.l I2OO6 ) the smoothing length h has been 
considered a function of time by relating it to the lo- 
cal particle density. A spatial and temporal smoothing 
length together with an appropriate symmetrization con- 
cerning particl e pairs hav e also been proposed ( Evrard 1988 



Hernguist k. K atz 198^ Nelson &: Papa,loizou _ 19931. 
Fulbright et al..,1995. : , Shapiro et al.ll 19961 : lOwen et al.l 



Liu et al-lbOO^ 



1994 



1998 



In original 3D ASPH hi varies in space and time. Sym- 
metry in both i, j indexes is widely adopted, where the evalu- 
ation of a symmetrized hij = {hi + hj)/2 and a symmetrized 
Kernel Wij = {Wij{hi) + Wji{hj))/2 are required according 
to: 



h 



n + l 



hi 



1/3 



(A14) 



where indexes n and n+l refer to time-step 
iHernauist fc Katzlll989l: iNelson fc Papaloizoulll99l | l994 : 



Fulbright et allll995l : IShapiro et al.lll996l : IOwen et al.lll998 : 
Liu et al] 20061 ). Such a choice is widely considered better 
than: 



/i7 = h- — 



1/3 



(A15) 



where h° and p° refer to initial values at time zero. Such 
a preference is due to the fact that because of non-linearity, 
instabilities can easily be produced es pecially in anisotropi c 
volume changes and fiow distortion ijMivama et al.lll984h . 
Equivalently, a further equation able to compute the "new" 
h at time-step n + 1 from the "old" h at time-step n 
llFulbright et al.lll995l : IShapiro et al]|l996l : lOwen et al.ll 19981 : 
iLiu et al.ll2006l ) is: 



I n+l I n 

h, ^ = hi 



or, by considering the continuity equation (1): 



h7^ 



h" 



1 + 



1 dpi 
Pi dt 



At" 



(A16) 



(A17) 



whose integration over time gives eq. (A14). This equa- 
tion is easily obtained by performing the derivative of the 



equation ph^ = const, expressing the conservation of parti- 
cle mass: 

dp/dt + 3ph^dh/dt = 0, dh/dt = -{h/3){dp/dt) / p, etc.. 
However eqs. (22), (23) or (2 4) ar e more convenient than 
eq. (A14). IShapiro et all (| 19961 ) and lOwen et~all (|l998h . pro- 
posed an adaptive method splitting the 3D scheme into 
three ID schemes formulating a factorized Gaussian Ker- 
nel of three ID Gaussian components. In such a scheme a 
tensorial computation of SPH equations has been developed 
and each ASPH particle enlarges or contracts as a spheroid 
rather than a spherule. They successfully applied their tech- 
nique to a shock front cosmological problem where ASPH 
spheroids give a better shock resolution compared to typ- 
ical SPH spherule without adopting any arti ficial viscosity 
term. In a further paper (|Owen et al.l Il998l ) the authors, 
admitting that artificial viscosity terms are necessary, es- 
pecially in the momentum equation, handle such artificial 
viscosity terms suppressing or turning on them according to 
some physical circumstances (mainly in rarefaction condi- 
tions). A technique tu rning on/off the artific i al vis cosity has 
also been described in lMorris fc MonaghanI (|l997h . 

ASPH models adopt the SPH same formulation, where 
either: 



'iij = ^{hi + hj) 
W,,=W(n,,h,,), 



(A18) 



instead of SPH h,, and W^j = W{r^j,h) l|Evrardlll988l ). 



Wij = \{W^j,^ + Wr, 
Wij,. = W{rij,h,), 



(A19) 



instead of SPH Wij — W{rij,h), are adopted 



(IHernauist fc Katz|[T989l ). The second formulation is mostly 

more currently adopted. 

Non-isotropic AS PH (|Shapiro et all 1 19961 : lOwen et all 
1 19981 : iLiu et all I2OO6I ) adopt an anisotropic algorithm to 

compute ellipsoid particle deformation and, consequently, 
the anisotropic smoothing length, according to the local par- 
ticle concentration. Such a scheme is mainly used in simu- 
lations of 2D and 3D oblique shocks and of contact fiuid 
surfaces. The algorithm computes the element h^^i, where 
a,l3 — X, y, z, of the 3x3 symmetric matrix: 



in 



1 + 



0.5 / dv. 



dViE 



At" 



(A20) 



where hafti = hpai, is the projection of the ellipsoid 
characteristic semiaxes on the cartesian axes. The eigenvec- 
tors of the matrix are the directions along the three axes of 
the ellipsoid and the corresponding eigenvalues are the di- 
mensions of the ellipsoid along each axis. The determinant 
of the same matrix determines the normalization volume of 
each particle. 

The SPH conversio n of eq. (A20), si r irilarl y to the SPH 
expression of the V • v (|Monaghanlll985l . Il992l ) is: 



7 n 

It-aBi 



, 0.5Ai" sr^rrij / „ 



i=i 



VpijViaWij 



(A21) 
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A2 Conservative ASPH formulation 



iNelson fc Papaloizoul j 1991 1 19941 ) showed that energy con- 
servation improves if d/dh are introduced into both SPH 
momentum and energy equations. The inclusion of such 
terms modify substantialfy those equations in a non 
practical form. The forma l difficulties were overcome by 
ISprineel fc Hernauisd l|2002l l who derived an effective ASPH 
conversion of the pressure gradient contribution in the mo- 
mentum equation (eq. 2), conserving energy and entropy, 
according to the conservative ASPH equation: 



fi l^vz 



(A25) 



where f2ij, includes artificial viscosity terms. In con- 
servative ASPH approach, it is easy to update the particle 
smoothing resolution length hi, fixing the number of parti- 
cle neighbours. In fact, according to the SPH interpolation 
criterion, particle concentration rii — 'Y^'^-i ^ij- We remind 



Vij/hij. Therefore, if Nndgh represents the fixed number of 
neighbours, = {Nneigh/niY^^ . 



(A22) 



APPENDIX B: TESTS 



where fi = ( 1 -f 



hj dpi ' 



hi) and 



Hij refers to the artificial viscosity contribution. Smoothing 
length h was computed requiring that a fixed mass is con- 
tained within a smoothing volume: (47r/3)/ii pi = Mij where 
Mi^j = rrijNi^j refers to the global mass of Nij neighbours 
related to the ith particle. Each particle neighbour has a 
rrij mass. No further modifications to the energy equation 
are required. In a further paper (iMonaghan 2002) similar 
conclusion, as far as both SPH and XSPH methods are con- 
cerned, were reached with the aim of achieving better energy 
and entropy conservation. 

The dpi/dhi term is easily connected to the dWij/dhi 
by the simple relation: 



In this section, results of some tests are here reported re- 
garding models where either ID shock problems, or 2D free 
edge, or 2D transport themes have to be taken into account 
to respect the argument declared in this paper. Compari- 
son among GASPHER, SPH and ASPH numerical results 
are reported as well as theoretical analytical ones, whenever 
the theoretical analytical solution is known. The particle 
smoothing resolution length h, normally adopted through- 
out, is ft = 5 • 10~^ (in ASPH as the initial value), but than 
when explicitly written. 7 = 5/3 throughout. Once stated 
the validity of GASPHER for shock coUisional modelling, a 
particular attention is addressed both to free edge and to 
radial transport results regarding the main argument of this 
paper. 



dpi 
dhi 



N 



dm, 

' dhi 



(A23) 



where the derivative dWij/dhi strictly involves also the 
derivative of the h~^ in 3D as: d{Wijh~''^) / dhi. In this 
scheme, the conservative ASPH conversion of the Navier- 
Stokes equation (eq. A7) is: 



m y p- p^ 



Pi 



(A24) 



As far as the conservative ASPH energy balance equa- 
tion for the total energy E is concerned, 



Bl ID Sod shock tube tests 

In this section a comparison of an alytical an d GASPHER 
ID inviscid shock tube test results l|Sodlll978l ). is made. No- 
tice that the so called analytical solution of the Id shock 
tube test is obtained through iterative procedures left-right, 
applying to the discontinuity the Rankine-Hugoniot "jump" 
solution. Figg. Bl and B2 display results concerning the par- 
ticle density, thermal energy per unit mass, pressure and ve- 
locity, after a considerable time evolution at time T = 100. 
The whole computational domain is built up with 2001 par- 
ticles from A = to A = 100, whose mass is different, 
according to the shock initial position. At time T = all 
particles are motionless. 7 = 5/3, while the ratios pi/p2 =3 
and ei/e2 = 2, and pi/p2 = 3 and ei/e2 = 1 as displayed at 
the edges of Figg. 4 and 5, between the two sides left-right. 
The first 5 and the last 5 particles of the ID computational 
domain, keep fixed positions and do not move. The choice of 
the final computational time is totally arbitrary, since the 
shock progresses in time, = at the beginning of each 
simulation. Hence, the adimensional temporal unity is cho- 
sen so that dx/cs — 1. Being the sound velocity initially 
constant, this mathematically means / = c^. SPH results, 
adopting the same initial and boundary conditions, as well 
as the same particle smoothing resolution length h, together 
with the analytical solutions are also displayed in the same 
plots. 

Our GASPHER results, are in a good comparison with 
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Figure Bl. ID shock tube tests as far as both analytical (solid line) and both SPH and GASPHER (short dashes) results are concerned. 
Density p, thermal energy e, pressure p and velocity v are plotted at time T = 100. The initial velocity is zero throughout. 
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Figure B2. ID shock tube tests as far as both analytical (solid line) and both SPH and GASPHER (short dashes) results are concerned. 
Density p, thermal energy e, pressure p and velocity v are plotted at time T = 100. The initial velocity is zero throughout. 



the analytical solution. Discrepancies involve only 4-^5 par- 
ticle smoothing resolution lengths at most. This means that, 
GASPHER interpolations are effective in the case of shock 
collision case in so far as the Mach number flows regard the 
weak shock regimes when the Mach number ranges within 
[0, 2] at the first instant. The decrease of the particle smooth- 
ing resolution length could improve the whole result in so 
far as the artificial viscosity term (depending onh- eq. A12) 
is able to prevent particle interpenetration. 



B2 ID Blast wave 

Whenever in a shocktube the ratios pi /pi = ei/e2 S> 1, and 
consequently pxj pi — 1, and vi — V2 = 0, such a discontinu- 
ity is called a "blast wave". Being w = at the beginning of 
each simulation, the adimensional temporal unity is chosen 
as previously written in the ID Sod shocktube test before. 
In such a situation, the Mach number spans from to > 1 
values up to 10 -i- 20 or more at the first instant. Fig. B3 
displays a comparison of SPH and GASPHER results with 
the so called analytical solution, after a considerable time 
evolution at time T = 5. The analytical solution is con- 
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Figure B3. ID blast wave tests as far as both analytical (solid line) and both SPH and GASPHER (short dashes) results are concerned. 
Density p, thermal energy e, pressure p and velocity v are plotted at time T = 5. The initial velocity is zero throughout. 



sidered corrected in so far as pi/p2 ^ (7 + l)/(7 — 1). In 
the blast wave test here considered, pi/pi — ei/e2 ~ lO**, 
while other spatial, initial and boundary conditions, as well 
as the particle spatial smoothing resolution length are iden- 
tical to those chosen in the previous test. Fig. B3 displays 
that SPH and GASPHER results globally compare with each 
other and that they also compare with the analytical solu- 
tion wherever pi/p2 ^ (7 -|- l)/(7 — 1), that is wherever the 
Rankine-Hugoniot jump conditions hold. Beyond this limit, 
even the so called analytical solution is considered incor- 
rect. Being 7 = 5/3, the comparison is meaningful within 
Pi/p2 4. GASPHER profiles suffer of a lesser instability 
in those regions where particle concentration is larger, close 
to discontinuity profiles, where horizontal plateaus are more 
regular. In particular such behaviour can be addressed to 
the absence of any particle pair instability because of the 
analytical expression of the GASPHER adopted Kernel and 
to its radial spatial derivative. 



B3 2D expansion of the free edge of a squared box 

The test here discussed does not have an analytical solu- 
tion. However, it is interesting because it shows how pressure 
forces push away the free edge of the fluid computational 
domain, without any explicit dissipation, according to the 
chosen interpolation Kernel. Being in permanent gradual ex- 
pansion, any artificial viscosity contribution is statistically 
turned off, apart some contribution due to the shear flow 
close to the two marginal vertical fixed edges. 

The box is a square 4.8 x 4.8, having three fixed sides: 
two vertical unlimited sides (left - right), at X £ [0.,0.05] 
and at X G [4.75, 4.80], and the horizontal one at the bot- 
tom from Y = 0. to y = 0.05, while the fourth side at 
Y = 4.80, is free to expand towards the outer space. Parti- 
cles, whose mass rrii = 3 ■ 10~® are regularly located so that 
their mutual separation equals h. The initial thermal energy 
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Figure B4. SPH, ASPH and GASPHER XY plots of the 2D 
expansion of a portion of the free edge on the top of a squared 
box. Time T is shown. 



is e = 1.9188-10"^, while the initial Vi = throughout. The 
three fixed edges are composed of two lines of fixed particles, 
whose velocity is abruptly put to zero time by time. Notice 
that the above mentioned constraints have to be considered 
as geometric conditions not pertinent only to the edge par- 
ticles of the square box. In this way, particles can move only 
toward the Y direction, from Y = 4.8 onwards, and any 
particle horizontal translation is mechanically prevented. 

Since an analytical solution is unknown, we pay atten- 
tion to the conservation of the total energy Etot = + e 
per unit mass averaged for each particle and, at the same 
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Figure B5. Plots of specific total energy Etot averaged for each 
particle for the SPH, ASPH and GASPHER simulations of the 
2D expansion of the free edge in a box. Etot values are expressed 
in 10~^ units. Time T is reported on a logarithmic scale. 



time to the regular face of the expansion of the free hori- 
zontal edge at the top. Fig. B4 displays the advance of the 
free front at three selected times for the SPH, ASPH and 
GASPHER simulations. SPH and particularly ASPH fronts 
are without any doubt more advanced than the GASPHER 
front. This effect is the result of an incorrect computation 
of the pressure forces on the free edge of the computational 
domain as discussed in §5. This conclusion is stressed not 
only by the fact that the GASPHER flovir is more regular 
and free from defects, but also, as it is shown in Fig. B5, 
by the fact that the total energy per unit mass Etot is much 
better conserved than in the other two cases. As an order 
of magnitude, the degradation of the total energy is ~ 10~* 
for a totality of ~ 10* particles after a time T ~ 2 • 10* for 
both SPh and ASPH. This involves that on a single particle, 
dE/dt ^ 5-10"^^. Instead, in GASPHER this energy degra- 
dation is ~ 10^ times smaller. This implies that the choice 
of the interpolation Kernel is crucial in the conservation of 
prime integrals. 
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Figure B6. XY plots of ring density contour maps. Times are 
reported for each configuration (SPH or ASPH or GASPHER). 




B4 2D radial spread and migration of a Keplerian 
annulus ring 

The 2D radial spread and migration of a n isothermal Keple - 
rian annulus ring is widely described in lFrank et al.l (|2002l ) 
in the case of a constant physical viscosity v. At time T = 0, 
the surface density, as a function of the radial distance r, is 
described by a Dirac 5 function: S(r, 0) — M5{r — ro) /2Tiro, 
where M is the mass of the entire ring and To is its initial 
radius. As a function also of time, the surface density is com- 
puted via standard methods as a function of the modified 
Bessel function Iin{z): 



^1/4^ 



-h/i{2x/T) 



(Bl) 



Figure B7. Surface density S(r, r) in 10~^^ units as a function 
of radial distance from initial configurations at r = 0.018. Subse- 
quent snapshots are reported for each configuration both SPH or 
ASPH or GASPHER. 



where x = r/ro, r = 12vT/rZ'^. J^'E{x,T)dS = 
27r J ^{x, T)dr = const equals the annulus mass throughout. 
Time is normalized so that T = 1 is the Keplerian period 
corresponding to the ring at ro = 1. Example s of SPH vis- 
cous sp r ead on this argument ca n be found in Flebbe et al.l 
(|l994bl ): ISpeith fc Riffei^ (Il999l ): ISpeith fc Kiev! (|2003D . as 
well as in Costa et al. I (|2009l ') in SPH physically inviscid 
hydrodynamics on the basis that the shear dissipation in 
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non viscous flows can be compared to physical dissipa tion 
jMolteni et al.lll99ll : lMurravlll996l : lOkazaki et"ai1l2002l '). In 

particular an exhaustive comparison can also be found in 
iLanzafamd ([2008b, 2009;). 

In a non viscous particle Lagrangian fluid dynamics, 
any deviation from the initial strictly Keplerian kinematics 
is incorrectly due to the activation of artiflcial viscosity dis- 
sipation in the shear flow when two particles approach each 
other. This is an unavoidable consequence of the fact that 
dissipation is currently used to handle the direct head-on 
collision between pair of particles. To establish whether the 
adopted Kernel has a signiflcant role in the spatial trans- 
port phenomena, any gas pressure force component must 
be removed leaving active only the artificial viscosity dissi- 
pation in the momentum equation in a strictly isothermal 
fluid dynamics. Thus, the whole flow should keep its Ke- 
plerian behaviour because pressure forces are artiflciously 
erased, in so far as artiflcial viscosity dissipation stays in- 
active. In this test, the two marginal edges of the annulus 
(the inner and the outer ones) are considered as FE bound- 
aries. Adopting the same artificial viscosity formulation and 
the same parameters (a = 1 and P = 2 - see App. A) both 
for SPH and for ASPH and for GASPHER simulations, our 
aim is to check which technique shows the smaller devia- 
tion from the initial Keplerian tight particle distribution in 
isothermal conditions, keeping constant both the sound ve- 
locity and the specific thermal energy. SPH-derived tech- 
niques turns on the artificial viscosity dissipation whenever 
two close particles approach with each other. This happens 
also for shear fiows. However for inviscid ideal shear flows 
this is an incorrect result without any gas compression. 

A signiflcant comparison of GASPHER to SPH and to 
ASPH is displayed in Fig. B6, where XY density contour 
map plots are shown at the same r. The radial distributions 
of surface density are displayed in Fig. B7, according to the 
restricted hypotheses of the standard mechanism of physical 
di ssipation (constant diss i pation, zero initial th ickness). As 
in lSoeith fc RiffertI (jl999l ): ISpeith fc Kiev! (|2003l ). the initial 
ring radius is at ro = 1, whose thickness is Ar = 0.5, is com- 
posed of 40000 equal mass (m^ — 2.5 • 10~^^) pressureless 
Keplerian (v = VKepi, V • v = at T = 0) SPH particles, 
with h = 9 - 10~^, with Cs = 5 - 10~^, and with initial density 
radial distribution corresponding to the analytical solution 
at time T, whose r — 0.018. To this purpose, a random num- 
ber generator has been used. The central accretor has mass 
normalized to M = 1. The k i nematic shea r dissipation is 
estirnated jMolteni et al.lll99ll : lMurravl[T99^ : lOkazaki et"al] 
l2002l : lLanzafamdl2008bl . l2009l ) as « Csh. 

GASPHER radial spread is without any doubt the nar- 
rower one, while the ASPH one is naturally the larger be- 
cause of the increasing particle smoothing resolution length 
h affecting the artiflcial viscosity analytical expression. For 
this reason, its physical dissipation counterpart u ^ Csh can- 
not be kept constant even preventing any disc heating. Only 
in the case of ASPH modelling, the initial h value is assumed 
to compute r. Notice that these results are obtained accord- 
ing to the correlation v x Csh ijMolteni et al.|[l99ll ) in the 
expression where r = 12i^T'/r^'^, which appears as t he most 
appropriate. In fact, cons idering u ^ O.laspHCsh (|Murravl 
1 19961 : lokazaki et al.l[20oi ) with asPH ~ 1, it is necessary 
a time T ten times longer to get the same t. This involves 
a larger annulus spread as far as the numerical results are 



concerned. According to this results, the Kernel choice is 
determinant also in the generation of kinematic deviations 
from the initial Keplerian distribution due to the incorrect 
SPH dissipation because of the particle shear approaching 
in the non viscous ideal flows. 

Notice that the density radial distribution, as far as the 
ASPH modelling is concerned, better flts the spread of the 
theoretical radial distribution (here not represented). This 
is a fair result in so far as we are interested in determining 
which artificial dissipation, coupled with the choice of the 
interpolation Kernel, determines a density radial profile, to 
be compared with the theoretical one, when a physical dis- 
sipation 1/ is considered. However, this is another aspect, 
regarding the study of either the physical dissipation in a 
viscous fluid dynamics or its artificial numerical dissipation 
counterpart in a non viscous approach, which is far from the 
scope of the test here proposed. 

Circular rings, appearing in Fig. B6 for both numerical 
schemes, are an unavoidable effect due t o the Lagrangian 
artic l e-based technique, as discussed in ISpeith fc RiffertI 



iq ue, as 

iiiooi) 



19991 ) : ISpeith fc Kle 



REFERENCES 

Belvedere, G., Lanzafame, G., Molteni, D. 1993, AfcA, 280, 
525 

Benz, W., Bowers, R.L., Cameron, A.G.W., Press, W. 1990, 
ApJ, 348, 647 

Bonet, J., Lok, T.S.L. 1999, Comp. Meth. in App. Mechan- 
ics and Engineering, 180, 97 

Bonet, J., Kulasegaram, S., Rodriguez-Paz, M.X., Profit, 
M. 2004, Comp. Meth. in App. Mechanics and Engineer- 
ing, 193, 1245 

Boris, J.P., Book, D.L. 1973, JCoPh, 11, 38 

Costa, v., Pirronello, V., Belvedere, G., Del Popolo, A., 
Molteni, D., Lanzafame, G., 2010, MNRAS, 401, 2388 

Drimmel, R. 1996, MNRAS, 282, 982 

Evrard, A.E. 1988, MNRAS, 235, 911 

Feldman, J., Bonet, J. 2007, Int. J. Num. Meth. Eng., 72, 
295 

Flebbe, O., Miinzel, H., Riffert, H., Herold, H. 1994a, Mem. 

S.A.It, 65, 1049 
Flebbe, O., Miinzel, H., Herold, H., Riffert, H., Ruder, H. 

1994b, ApJ, 431, 754 
Frank, J., King, A.R., Raine, D.J., 2002, "Accretion power 

in astrophysics", Cambridge Univ. 
Fulbright, M.S., Benz, W., Davies, M.B. 1995, ApJ, 440, 

254 

Hernquist, L., Katz, N. 1989, ApJS, 70, 419 
Kaisig, M. 1989, AfcA, 280, 525 

Katz, N., Weinberg, D.H., Hernquist, L. 1996, ApJSS, 105, 
19 

Imaeda, Y., Inutsuka, S 2002, ApJ, 569, 501 
Lanzafame, G. 2003, AfcA, 403, 593 

Lanzafame, G. 2008a, "The Role of Physical Viscosity in 
Accretion Disc Dynamics in Close Binaries and AGN", 
"Numerical Modeling of Space Plasma Flows/ Astronum 
2007", N. V. Pogorelov, E. Audit, and G. P. Zank eds., 
ASP Conf. Series 385, p.ll5 

Lanzafame, G. 2008b, PASJ, 60, 259 

Lanzafame, G. 2009, AN, 330, 843 



18 G. Lanzafame 



Lanzafame, G., Belvedere G., Molteni D. 1992, MNRAS, 
258, 152 

Lanzafame, G., Belvedere G., Molteni D. 1993, MNRAS, 
263, 839 

Lanzafame, G., Maravigna F., Belvedere G. 2000, PAS J, 

52, 515 

Lanzafame, G., Maravigna F., Belvedere G. 2001, PASJ, 

53, 139 

Lanzafame, G., Belvedere G., Molteni, D. 2006, A&A, 453, 
1027 

Lasota, J.R 2001, New Astr. Rev., 45, 449 

Lattanzio, J.C., Monaghan J. J., Pongracic, H., Schwarz, 

M.P., 1985, MNRAS, 215, 125 
Liu, M.B., Liu, G.R., Lam, K.Y. 2006, Shock Waves, 15, 

21 

Lubow, S.H., Shu, F.H., 1975, MNRAS, 198, 383 
Meglicki, Z., Wickramasinghe, D., Bicknell, G.V. 1993, 

MNRAS, 264, 691 
Miyama, S.M., Hayashi, C., Narita, S. 1984, ApJ, 279, 621 
Molteni, D., Belvedere, G., Lanzafame, G. 1991, MNRAS, 

249, 748 

Monaghan, J.J. 1985, Comp. Phys. Rept., 3, 71 

Monaghan, J.J. 1992, ARA&A, 30, 543 

Monaghan, J.J. 1994, JCoPh, 110, 399 

Monaghan, J.J. 1997, JCoPh, 136, 298 

Monaghan, J.J. 2002, MNRAS, 335, 843 

Monaghan, J.J. 2006, MNRAS, 365, 199 

Monaghan, J. J., Lattanzio, J.C. 1985, A&A, 149, 135 

Morris, J. P., Monaghan, J.J. 1997, JCoPh, 136, 41 

Murray, J.R. 1996, MNRAS, 279, 402 

Nelson, R.P., Papaloizou, J.C.B. 1993, MNRAS 270, 1 

Nelson, R.P., Papaloizou, J.C.B. 1994, MNRAS 265, 905 

Okazaki, A.T., Bate, M.R., Ogilvie, G.L, Pringle, J.E. 2002, 

MNRAS 337, 967 
Owen, J.M., Villumsen, J.V., Shapiro, P.R., Martel, H. 

1998, ApJSS 116, 155 
Pringle, J.E., Verbmit, F., Wade, R.A. 1986, MNRAS, 221, 

169 

Savonije, G.J., Papaloizou, J.C.B., Lin, D.N.C. 1994, MN- 
RAS, 268, 13 
Sawada, K., Matsuda, T. 1992, MNRAS, 255, 17 
Sawada, K., Matsuda, T., Inoue, M., Hachisu, I. 1987, MN- 
RAS, 224, 307 
Shakura, N.l. 1972, Astron. Zh., 49, 921 
Shakura, N.I. 1973, SvA, 16, 756, Engl. Transl. 
Shakura, N.L, Sunyaev, R.A. 1973, A&A, 24, 337 
Shapiro, P.R., Martel, H., Villumsen, J.V., Owen, J.M. 

1996, ApJSS, 103, 269 
Sod, G.A. 1978, JCoPh, 27, 1 
Speith, R., Kley, W., 2003, A&A, 399, 395 
Speith, R., Riffert, H., 1999, JCoAM, 109, 231 
Springel, V., Hernquist, L. 2002, MNRAS 333, 649 
Spruit, H.C., Matsuda, T., Inoue, M., Sawada, K. 1987, 

MNRAS, 229, 517 
Vaughan, G.L., Healy, T.R., Bryan, K.R., Sneyd, A.D., 
Gorman, R.M. 2008, Int. J. Num. Meth. Fl., 56, 37 



